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We investigate finite number effects in collisions between two states of an initially well known number of 
identical bosons with contact interactions, oscillating in the presence of harmonic confinement in one dimen- 
sion. We investigate two N 12 (interacting) ground states, which are initially displaced from the trap center, and 
the effects of varying interaction strength. The numerics focus on the simplest case of A' = 4. In the non- 
interacting case, such a system would display periodic oscillation with a half harmonic oscillator period (due to 
the left-right symmetry). With the addition contact interactions between the bosons, collisions generate entan- 
glement between each of the states and distribute energy into other modes of the oscillator. We study the system 
numerically via an exact diagonalization of the Hamiltonian with a finite basis, investigating left/right number 
uncertainty as our primary measure of entanglement. Additionally we study the time-evolution and equilibra- 
tion of the single-body von Neumann entropy for both the attractive and repulsive cases. We identify parameter 
regimes for which attractive interactions create qualitatively different behavior to repulsive interactions, due to 
the presence of bound states (quantum solitons) and explain the processes behind this. 
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I. INTRODUCTION 

Dilute gases of Alkali atoms have proved a powerful tool 
for the experimental investigation of quantum mechanical 
phenomena, from the level of single atom physics up to to 
mesoscopic levels via the creation of Bose-Einstein conden- 
sates (BECs) [1,2]. Much of the interest stems from the abil- 
ity to experimentally realize many theoretically interesting po- 
tentials, such as optical lattices [3], double well potentials [4] 
and periodic kicking [5], with the ability to control the effec- 
tive dimensionality and interaction strength via Feshbach res- 
onances. Another interesting property is the ability to support 
both bright and dark solitons [6, 7]. 

Experimentally it is possibly to tune i-wave scattering 
lengths to both positive and negative values [8-10]. How- 
ever, beyond a certain critical number (which is dependent 
on the trapping configuration and scattering length), the nega- 
tive scattering length (attractively interacting) systems are un- 
stable to collapse [11-14]. If trapping potentials are present 
in two spatial dimensions, attractive condensates can exhibit 
self trapping, i.e., localization (at least in terms of pair cor- 
relation functions) in a direction free of external potentials. 
In quasi one-dimensional (ID) geometries, attractive BEC's 
form bright matter-wave solitons with particle-like dynamics 
for the center of mass [15]. Parameter regimes for which 
systems are quasi- ID have been investigated via variational 
techniques [13], along with effective potential approximations 
to deal with residual 3D effects [16], leading to higher order 
effective non-linearities. In addition to this, bright gap soli- 
tons [17] have been created from repulsive atoms in optical 
lattices by exploiting anomalous dispersion to give the atoms 
a negative effective mass. 
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Negative scattering lengths also give interesting possibili- 
ties in double-well and lattice physics. Repulsive interactions 
between atoms are known to give rise to the famous Mott in- 
sulator state [18], with a near definite atom number per lattice 
site. If one has a definite number of atoms per site, there is 
effectively a total uncertainty in relative phase between lattice 
sites and thus no phase coherence. A measurement of relative 
phase should give totally random results and indeed this is 
what one finds when imaging the moment distribution of such 
a lattice, no distinguishable interference patterns. Attractive 
interactions could in theory be used to squeeze number statis- 
tics the opposite way, such that the ground state would tend 
to a superposition of a quantum soliton (A^ atom bound state) 
delocalized over every lattice site. When only two sites are 
present, such a state is referred to a NOON state [19], which 
is useful for non shot noise limited interferometry [20]. How- 
ever, systems where the ground state is such a superposition 
are known to be extremely unstable to temperature, as phase 
differences between the two sites have almost no energy cost, 
thus typically replacing quantum uncertainty with statistical 
uncertainty. It is therefore preferable to create such states dy- 
namically, for example by splitting a moving quantum soli- 
ton [21, 22]. 

Any closed quantum system with no decoherence effects 
will be described by a wavefunction that will evolve determin- 
istically. As such the wavefunction at any point in time 
maps back to a unique |i/'(0)). Recent experiments have shown 
great possibility to observe this deterministic behavior in sys- 
tems with a small number of cold atoms [23, 24], with dynam- 
ics that can be analytically calculated and with precise tun- 
ing available in the scattering length and confinement poten- 
tials. Strongly correlated effects and quantum superpositions 
are generally much easier to achieve in few-body systems. 
Despite this one can still envisage collective properties (such 
as expectation values of operators) of a time-dependent finite 
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system tending to constant values when averaged over rea- 
sonable timescales, or relaxation of local operators, as shown 
in [25]. Non-integrable systems, upon coupling to another 
larger system, usually tend to an equilibrium configuration at 
long times, independent of the initial state of either system 
(except for the total energy); however recent theoretical ob- 
servations have thrown doubt on this [26]. Additionally, when 
two coupled systems contain a similar number of elements 
the situation is less clear still. Our system is non-integrable 
and contains two initially independent subsystems of the same 
size; hence, we are interested to what extent equilibration oc- 
curs or where it is resisted. Quantum systems, for example 
atoms populating sites in an optical lattice [27] are known to 
show partial revivals of the initial state in time, but are gener- 
ally observed to show weaker revivals as time progresses in an 
apparent damping; we are interested in whether certain mea- 
sures, specifically the number to the left and right of are trap 
and the single body von Neumann entropy, tend to constant 
values when averaged over sufficient timescales. 

This paper is organized as follows: Section II introduces 
the one dimensional Hamiltonian and the unit rescaling to har- 
monic oscillator lengths, constant throughout the paper Next, 
the initial condition is introduced, with specific cases of inter- 
est mentioned. Section III discusses observables and measures 
of entanglement that we will use to investigate the system, in- 
cluded the variation in the number to either side of the trap 
center and the single-body von Neumann entropy. Section IV 
begins an analytic investigation of the system, focusing on the 
mechanisms by which interactions modify the dynamics of 
each displaced state and generate entanglement. Section V 
discusses a possible experimental realization of the system, 
using ultracold atoms in an optical lattice, with parameters 
discussed for Cesium. Section VI contains a brief descrip- 
tion of the numerical method, based on exact diagonalization. 
Section VII presents numerically obtained results for the evo- 
lution of our observables and entanglement measures in the 
system. Section VIII summarizes and concludes. 



where M is the mass and the axial (angular) trapping fre- 
quency; assuming a radial trapping frequency of oj^, the cou- 
pling parameter gi^, - Ihcj^Os, with a, the (3D) i-wave scat- 
tering length [31, 32]. A satisfactory condition for this Hamil- 
tonian to be valid is N\as\ y/h/Mcji- and ksT <g: hcL),-, how- 
ever it is likely to be still be valid for ksT ~ ticor, i.e., as long 
as thermal excitations are unlikely to significantly populate 
radial modes. 

We use harmonic oscillator units (codified ?& fi - - 
M = I), meaning that length is in units of ^fijMujx, time in 
units of I/w.v, and energy in units of hoj^; a harmonic oscilla- 
tor period is then 2-k. The Hamiltonian rescales to 



^{x), (2) 



where g — gio -^^mJWoj^ is the new dimensionless coupling 
parameter' In first quantization we can express this same 



Hamiltonian (for particles) as 



Idxl 



N k-\ 

+ ^ ^ ^ S{xk - Xj) , (3) 

k=2 j=l 



where x^ are the coordinates of the individual particles (gen- 
erally considered to be ultracold atoms), and is a shorthand 
for the set of all coordinates {xi,X2, . . . , x^]. As the external 
potential is harmonic, H{x) can be partitioned into two mutu- 
ally commuting components [28], one describing the center 
of mass (giving rise to the Kohn mode [33]), and the other 
describing the remaining degrees of freedom. This separa- 
tion can be exploited computationally, as the center-of-mass 
dynamics are those of a simple harmonic oscillator and can 
therefore be described exactly, reducing the effective dimen- 
sionality of the computational problem to N - I. 



B. Initial condition 



II. SYSTEM 



1. General N-body case 



A. Hamiltonian and unit rescaling 

We consider an effective one-dimensional (ID) system 
[taken to be reduced from a three-dimensional (3D) con- 
figuration where the radial degrees of freedom are strongly 
confined by a harmonic trapping potential] of structureless 
bosons subject to attractive or repulsive contact interactions 
y(|jt:i - X2I) = giD^ixi - X2), i.e., a Lieb-Liniger-(McGuire) 
gas [28-30], with the addition of an axial harmonic confining 
potential. In second-quantized form, this can be described by 
the following Hamiltonian: 



dx^'Hx)\ ^ + —\mx) 



dx^'HxyvHxyi'ixyi'ix) , (i) 



We consider a highly non-mean-field-like initial condition, 
taking two A^/2-atom ground states (for a given g), equally 
and oppositely displaced from the trap center by a distance 
xo, and symmetrizing. The initial (f = 0) wavefunction is then 



HX, 0) =^ J] f'^-Xxi - Xo, ..,XN/2 - xo) 



X f^'^KxN/l+l +Xo,...,Xn+Xo), 



(4) 



where /'"^^^(xi, . . . ,xn/2) is the ground state for N/2 atoms 
(generally numerically determined) in the harmonic trap, {P} 
is the set of all permutations of and B is a normalizing fac- 
tor Such an initial condition may be motivated by the idea 



This relates to the parameter y of [28] through y = [g(N - 1)] 
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of making two separate BECs and allowing them to collide 
within a harmonic trapping potential, or from rapidly modi- 
fying a Mott insulator state in an optical lattice (as we will 
discuss in Section V). If the left and right components are 
well separated, i.e., the width of the atomic density distribu- 
tion corresponding to is significantly less than x^, then 
there is a well defined number of N/2 atoms either side of the 
trap, and left- and right-atoms are distinct by virtue of their 
position. Furthermore, as the center-of-mass dynamics are de- 
coupled [28] and straightforward to determine, the dynamics 
experienced by an initial condition such as ip can be readily 
extended to incorporate any initial condition for the center of 
mass, e.g., in particular, an overall oscillation about the trap 
center [15]. 

Conveniently, t//{x, t) is in the ground state of the center- 
of-mass component of H(x). To show this, we first define 
(unnormalized) Jacobi coordinates, consisting of the center- 
of-mass coordinate 



1 



and - 1 further independent coordinates 



(5) 



1 

7=1 



(6) 



where k e {2, 3, . . . , A^). Using the Jacobi coordinates for 
N/2 particles, we can partition the A^/2-particle ground state 
into center-of-mass-dependent and independent components: 



^ ^ = fi^i, . . . , iN/2) e-^-'c<«/2)'r Substituting 
in Eq. (Al), we can then define /f^^^' through 



./4 



^N12 2/1 



(7) 



where /('^Z^' (as it can also be written as a function of 
{^2,^i,---,^N/2} only) is clearly independent of xc(n/2)- 
Within an expanded set of coor dinates, /'^^^j^^j^ _ 
is also clearly independent of xc(n), as is (by symmetry) 
f^'^l^\xNi2+\, . . . , xm). Noting further that displacement by xq 
will not affect that part of independent of the center-of- 
mass coordinate, then for the identity permutation of i/i 

f'''^'^\xi -Xq,..., Xn/2 - Xo)f'-'^'^\xN/2+l +Xo,...,Xn + Xq) 

^f^''l^\xu...,xr,i2)F^^\xNi2.u...,xr,) 



X e 



(8) 



By the identities Eq. (Al) and Eq. (A7), the exponential re- 
duces to e"'^'cw)^^e"'^(^" f'/'^--^o/2)^ j e., a term proportional 
to the center of mass ground state multiplied by a function 
of independent Jacobi coordinates. The identity permutation 
of i/r can thus be written as a product of the center of mass 
ground state and a function of the other independent coordi- 
nates. This separation off of the center of mass ground state 
occurs for every permutation of the coordinates Xk, and so we 



conclude that the center-of-mass component of t/r is indeed in 
the ground state. 

Taking a slightly different initial condition, when one com- 
bines ground states from two trapping potentials which are not 
equal to the final potential (with, e.g., tighter harmonic trap- 
ping), will introduce a breathing motion, which can still be 
considered separately from the remaining dynamics. It is also 
significant to note that the kind of initial condition we con- 
sider does not have a well defined relative phase between the 
left and right components [34]. If a relative number uncer- 
tainty between left and right were to develop then this would 
no longer be the case, and a meaningful relative phase could 
in principle be extracted. 



2. Time evolution for the non-interacting case 

If we take the case where g = 0, we can express the full 
time dependent wavefunction [which we label i^o(x, f)] ana- 
lytically, as a symmetrizing product of two A^/2-atom product 
states 

^ NI2 N 

i//Q{x, t) = — ^ ^ Y\ 'f'(^k, -Xu, f) Y\ 'f'(^j' -^0' ^) ■ (9) 

IP) k=l ;=«/9+i 



j=NI2+\ 



Here (p{x, +xq, 0) is a Gaussian displaced by ±xo from the trap 
center, and [15, 35] 



IW'^ I [x-xocos(f)]2 
(/)(x, xo, = I - I exp I 

X exp(![f/2 - Xq cos(f)x + Xq sin(2f)/4]) , 

corresponding to an energy-per-particle of £ = (x^ + Y)/ 
where the normalization constant Bq - \ + (9(e"^^o). 



(10) 



3. = 4 special case 

If N - A, the Z*^^* appearing in Eq. (4) are known analyti- 
cally [36, 37], and if g < may, for sufficiently large g and xq, 
be considered to be bound-state dimers, held within an overall 
harmonic trapping potential. The general form is given by 



/2)(xi,X2)-^f/l-v,l/2 



[Xi -X2]2\^_^2/2g-.4/2 



(11) 



= 7Vf/(-v, 1/2,^^/2) e-*5/^e 



with U Tricomi's confluent hypergeometric function, N a 
normalization constant, and v the effective quantum number 
(equal to zero for g - 0), as determined by the transcenden- 
tal equation r(l/2 - v)/T{—v) - -g/T?^'^. This state has an 
energy of 2v -H 1, where there is a contribution of 1/2 due to 
the center of mass. Equation (11) can then be inserted into the 
initial condition 

(/'(xi,X2,X3,X4,0) y /^'(Xi -Xo,X2-Xo) 

V4! ^ (12) 

X /^\X3 H- Xo, X4 H- Xo) , 
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where {!P) is the set of all 4! permutations of {jci, X2, X3, X4). 
Note that, as f^^\xi,X2) = f^^\x2, xi), the number of distinct 
permutations actually reduces to 4!/2!2! =6. 



III. OBSERVABLES AND MEASURES OF 
ENTANGLEMENT 

A. Left/right number 

For our system, one useful measure to track the generation 
of entanglement is the variance in particle number to the left 
and right of the system's center of mass (which we will gen- 
erally consider to be fixed at the origin). The initial condi- 
tion we consider has N/2 atoms to either side with essentially 
no possibility of, say, N/2 + 1 to the right and N/2 - 1 to 
the left (probabilities for measuring such unequal partition- 
ings decrease Gaussianly with the initial separation). Hence, 
the left- and right-particle-number-variance will initially be 
zero. As the left- and right-particles approach and collide, all 
number partitionings become possible, and so this measure is 
only informative when the particle density at the location of 
the center of mass is small. 

We define a "number-to-the-right" operator 

Nr = I dx4'\xmx) , (13) 
Jo 

[or in first quantization Yjk^i ®(xk), where is the Heaviside 
step function] ; imaging one side of the trap would correspond 
to a projective measurement into the eigenstates of this oper- 
ator, as is discussed in section V. The expectation value of Nr 
is the mean number of particles on the right-hand-side — as 
the system is parity preserving, (Nr) = N/2 for all time for 
the initial conditions we consider 

The more informative number-to-the-right variance is 

AA^« = <^^) - {NRf , (14) 

which, for our initial condition of two well separated left-and- 
right components of definite number, should be « 0. From 
Eq. (B6), the variance for a product state i/f(x) = Ylk=i 'P(^k) 
[symmetric about the trap center so that (Nr) - N/2] is 

ApNr = <^r)(1 - {Nr)/N) = N/4- , (15) 

which evaluates to unity if = 4 (this is however the same as 
a symmetric superposition of one and three atoms to the right 
/left). It can also be shown (Appendix B 1) that for the case 
of = 4 and no interactions {g = 0) [given by Eq. (9)], this 
variance evolves as 

ArNr = 1 - evfixo cos(f)) + 0(e-^'«) , (16) 

with erf the error function.' Hence, we have a function with 
period T - n, which is equal to unity when t - {n + l/2)n and 
vanishingly small in jcq when t - nn. 



In general our wavefunction is not an eigenstate of Nr, and 
contains components of different Nr eigenstates (for some 
given overall A^, meaning that an additional specification of 
number-to-the-left operator eigenstates is not necessary). 

One can however calculate expectation values of operators 
defined over restricted regions of state space, specific to hav- 
ing exactly n (of A^) atoms to the right of the trap center An 
expectation value for an operator O defined in this region is 
then 

dxi . . . dx„ J^^dx„+i . . . dxN t//*(x) 0{x)i//{x) 

{0)n,N-n - — -^j ■ 

Jg dxi... dxn dxn+\ ■ ■ ■ dxN \^{x)Y- 

(17) 

This is equivalent to taking the usual expectation value over a 
new wavefunction il/n,N-n{x) defined by 

^ n N 

l/Jn,N-n(^ = '^(-^Zn®^-^'^^ n ' ^^^^ 

yjPn.N-n p k=l >=n+l 

where !P is the set of all unique permutations, of which there 
are N\/n\{N - «)!, and P„^n-„ is a normalizing factor, giving 
the probability of finding nof N atoms to the right (or equiva- 
lently N-nio the left) of the trap center Each such wavefunc- 
tion is an eigenstate of .^Vr, with eigenvalue n. In principle one 
can partition the Hilbert space in such a way that it is the ten- 
sor product of a subspace describing only how many particles 
are to the left/right of the trap center, and a subspace describ- 
ing all other relevant properties of the system state. We may 
denote the set of eigenstates of Nr spanning this "number sub- 
space" by {|A^ - n,n}}, such that 

Nr\N -n,n) ^n\N -n,n) . (19) 

We can also study expectation values of a distance-to-the- 
right operator J^^ dx xWHxy^{x) [Yjk=\ ®ixk)xk in first quanti- 
zation] and associated higher-order moments, which will trace 
particle-like tracks (and widths around them) for state compo- 
nents of different right-hand number n. 

B. von Neumann entropy and relaxation 

Averaging over all individual particles results in the single- 
body density matrix 

pix, x, f) = {^'Hx'y^ix)) , (20) 

which is normalized to the total particle number A^. From this, 
single-body properties of the many-body system may be de- 
termined, specifically the von Neumann entropy^^ 

5vN(f) = -Tr{(p/A^)ln(p/A^)) . (21) 

Relaxation, in the sense of tending to states of higher en- 
tropy, is not present if the system is fully integrable, i.e., when 



^ Satisfying erf(0) = and erf(±x) — > ±[1 - exp(-A:^)/( ^x)} as jc — * 00. 



^ This is sometimes referred to as tlie Invariant Correlation Entropy 
(ICE) [38]. 
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g = 0, or if the trapping is removed and the eigenstates are 
given by the Bethe ansatz [30]. However, as the integrabihty 
is broken by the trapping, we expect some degree of thermal- 
ization due to (previously forbidden) mixing between states. It 
is of interest to determine how such thermahzation timescales 
vary with the interaction strength and initial separations. 

For a product state, p has a single non-zero eigenvalue of 
value A^, meaning S'vn — > (this is equivalent to a Bose- 
Einstein condensate being exactly described by a Gross- 
Pitaevskii wavefunction). A larger value of 5vn indicates oc- 
cupancy of multiple eigenstates of p, equivalent to population 
of non-condensate modes due to to thermal excitations, or to 
quantum or dynamical depletion [39, 40]. 

If the system equilibrates, 5 vn will tend to a constant value. 
As our initial conditions result in repeated collisions at the trap 
center, the value of S'vn shows distinct oscillations that decay 
only slowly. We therefore also consider a time average over 
an oscillator period 



S VN(f) 



2 pt+ljr 



dt' Ssnit') , 



along with its variance 



ASvnCO 



df 



S VN(f') 

2n 



■ S vn(0 



(22) 



(23) 



If 5vn(0 tends to a constant value, this will be shown by a 
relaxation of 5vn(0 to a constant value, and a relaxation of 
A5vn(0 to 0, with the relaxation of 5vn(0 tending to occur 
on a significantly faster time scale than that of A5vn(0- 



IV. ANALYSIS OF THE INTERACTING SYSTEM 



occurs infinitely more often than xi = X3, which necessar- 
ily implies X2 - X3 and so is a set of lower dimensionality. 
As [Hi,Hk] = 0, if we neglect Hj our system can be de- 
scribed by a tensor product of the left and right components.'* 
Each Hamiltonian Hl/r can further be split into center-of- 
mass ^^^^ and relative parts generating the dynamics 

of the left and right center-of-mass and relative coordinates 
[xc(L) = (xi + X2)/2, xc(R) = (x3 + X4)/2, xr(z.) = x2 - xi, and 
■^R(«) = X4-X3, respectively], which again mutually commute. 

We consider the center-of-mass wavefunction of an n atom 
cluster, which is a Gaussian displaced from the trap center by 
some value X„. Without the influence of Hj our system con- 
sists of two indistinguishable clusters (with internal degrees of 
freedom considered to be in the ground state) undergoing sim- 
ple harmonic motion. The primary reason for separating the 
Hamiltonian in this way is that our initial condition is in the 
ground state of H^^jl, and is a displaced ground state of 
hence any change to these wavefunctions is an excitation of 
the system. 



B. Perturbative introduction of Hi 

1. Overview 

We consider the effect of introducing the Hamiltonian Hj, 
from Eq. (24), to the system. We look at three notable effects: 
changes to the wavefunction describing the left/right separa- 
tion of the clusters; changes to the internal degrees of freedom 
within the clusters to the left and right; and interactions trans- 
ferring atoms from one side to the other, creating a symmetric 
superposition. 



A. Left-riglit separation of tlie Hamiltonian 



2. Inter-cluster wavefunction changes and pseudo-periodicity 



As our initial condition consists of left and right compo- 
nents which are well separated and therefore distinguishable, 
we can initially treat the left and right components separately. 
As these left and right clusters only interact for a short-time 
during collisions in the center (so long as they stay as distinct 
clusters), it makes sense to treat interactions between these 
clusters perturbatively at early times. We therefore split the 
Hamiltonian into three, restricting the coordinates to the re- 
gion xi < X2 < X3 < X4, which is sufficient due to Bose 
symmetry. The three components are 



The center-of-mass wavefunctions of each side, described 



'2 8x1 



Hi{x\_,X2) = ^ 

i/«(x3,x4) = 2. + y 

Hi(x2, X3) = g6(xi - xi) . 



+ g5{x2-x{) 



+ g6(M - x-i) 



(24) 



The reason only adjacent interaction terms \5{xk - xj) with 
k - j - 1] remain is that the other terms constitute a set of 
zero measure in the region we are considering, i.e., x\ - X2 



by Hf' 



H 



(C) 



g , can change, so long as the global center-of- 
mass wavefunction remains constant. Such changes lead to 
entanglement between the left and right clusters, to see this we 
note initially the two cluster wavefunction could be written as 
a product of left and right sides 



4ioixc(i),xc(R)) e 



-I- perm , 



(25) 



with "perm" denoting the permutation of R and L. This can 
be written in such a way as to explicitly separate the global 
center of mass: 

UxciL), xc(R)) = e-I*'"^^^'"!'^^ (26) 



* Commuting Hamiltonians imply exp(-![//i + Hi(\t)\il/) = 
exp(-iHLt)\i//[^) exp{-iHfit)\if/K), i.e., tlie time evolution operator can 
be separated. 
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The first term describes the global center-of-mass and is there- 
fore fixed, the latter term, however, will be modified by inter- 
actions. Any such change (other than modifying xq or multi- 
plying by exp(ip[xc(L) - xc(r)]), which are simply rescalings 
of the initial position and kinetic energy, respectively) means 
there will be terms involving products of the form xc(l)Xc(R), 
such that the wavefunction cannot be separated, indicating en- 
tanglement between the left and right sides. Such entangle- 
ment is notable in the context of solitons in free space, as in- 
tegrability means collisions cannot create entanglement once 
the states are asymptotically separated, although higher order 
non-linearities can also lead to entanglement [41]. Addition- 
ally, during collisions with attractive (repulsive) interactions, 
each cluster will accelerate (decelerate), subsequently return- 
ing to near its initial velocity, leading to a pseudo-periodicity. 



3. Intra-cluster wavefunction changes 



The internal degrees of freedom described by are ini- 
tially in the ground state. Interactions during collisions will 
introduce excitations, with the energy transferring from the 
center-of-mass energy of each cluster. By conservation of en- 
ergy this must reduce the amplitude of the oscillation. Attrac- 
tive interactions will suppress such excitations, as the energy 
separation between ground and first (even parity) exited state 
is greater than the harmonic oscillator level spacing, whereas 
for repulsive interactions this gap will be smaller. Note that 
when highly excited modes of the relative degrees of freedom 
xr(L),x^{R) are populated, these will always have a significant 
occupation for both L and R. One expects a qualitative differ- 
ence in behavior between the attractive and repulsive cases to 
occur when the change between the first and second relative 
excited states differs by an amount of order unity in harmonic 
oscillator units [ftw.v]- We note that for strongly attractive in- 
teractions,^, when Xq < -g/4 there is not enough energy to 
break the bound-state clusters, making the relative degrees of 
freedom effectively inaccessible, but this is beyond the scope 
of the present paper. 



4. Left/right atom transfer 

Finally, the interactions can transfer an atom from one side 
to the other, mixing to a set of states with a symmetric super- 
position of three and one atoms at either side of the trap (and, 
ultimately, back from this to the original state). There can- 
not be significant transfer to a state where there is a cluster of 
four atoms in the ground state (apart from the center-of-mass 
degree of freedom) on one side and zero on the side, due to 
the invariance of the center-of-mass wavefunction, unless the 
state has all four atoms directly at the trap center. The state 
satisfying this condition is the ground state of the system, and 
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FIG. I. (Color online): Energy difference Aiiij, = £3.1 -£2.2 between 
two+two and three+one atom ground state clusters, as a function 
of interaction strength in harmonic units. Analytic estimates from 
Eq. (27) are shown for comparison with the Tonks gas being the g — » 
00 limit. 



so the only possible population is that there at f = 0. Note 
that excited states of this four atom cluster do make up parts 
of the oscillating cluster states, it is simply a different basis to 
consider the system in terms of. 

A feature which distinguishes this effect from intra-cluster 
excitations is the energy difference between the two configu- 
rations, denoted AEint = £3,1 -£2,2- Forg < the ground state 
of a three atom relative Hamiltonian (that part of the Hamil- 
tonian independent of the center of mass) plus a single free 
atom is lower in energy than two sets of two atoms in their 
relative ground states. The opposite is true for g > 0, but the 
energy difference can only be of the order of the harmonic os- 
cillator energy spacings, and so suppression is unlikely unless 
xo is small. The energy difference AEint can take a variety of 
values when intra-cluster states are excited, but in the interest 
of studying transfer interactions, we look at the energy dif- 
ference between two isolated ground states of N = 2 atoms 
and an = 3 and N = I atom ground states. This can be 
estimated analytically in three limits: 



A£i, 



1 



if|gl« 1, 
if^ » 1, 
ifg « -1; 



(27) 



' In this regime energies scale as 
state [29]. 



l)/2 for an n atom ground 



the approximations used being overlapping non-interacting 
ground states, effective fermionization [42] (Tonks gas) and 
bound state clusters [29] with the first order energy correc- 
tion from the trapping potential [28], respectively. Numer- 
ically determined values of AEint are shown in Fig. 1; this 
energy proves to be an important quantity in the next section 
(note that this does not include the energy from the momen- 
tum/displacement of the clusters). Viewed classically, this 
transfer interaction causes transfer to a state where the kinetic 
energy of the clusters was different from the original by an 
amount equal to AEi„t, in order to conserve energy. 
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C. Mixing between diiferent number configurations via 
time-dependent perturbation tlieory 

We now investigate the atom transfer effect outlined in sec- 
tion IV B 4, predicted to be most significant for g <0. We can 
write our wavefanction at any point in time as 

\m) = C2,2(f)lfe(f)) + Ci,3(f)|lAl,3(0) + C0.4(0l'A(),4(f)) , (28) 

with \4'n,N-n(f)) normaUzed wavefunctions which are super- 
positions of states with n and N - n atoms to the left and 
vice versa, and {c„_Ar_„) a set of complex constants, the mod- 
ulus squares of which are the probabilities to find n or: N - n 
atoms either side. In order to qualitatively predict the incre- 
mental changes to {c„jv_„(f)) from before to after a collision, 
we use time dependent perturbation theory, assuming |g < 1 is 
a small parameter and neglecting any contribution from co,4(0 
(specifically at the time of collisions). We further assume the 
center-of-mass motion of each n,N-n atom cluster in |i^3.i(f)) 
undergoes harmonic oscillation and is periodic in time with 
period T - n, and that any internal relative excitations in both 
Vl'n.N-nit)) are small compared to the ground state. This ap- 
proximation is expected to work better for ^ < 0, for reasons 
outlined in Sec. IV B 3, and at short-times. As we initially 
have only C2,2 0, we assume |c3,i(f)l ^ k2,2(0l as a regime 
of validity. 

Formally, we perturb {Hi + Hr) by Hi [see Eq. (24)]. Our 
wavefunction 



W)) ^ ^,2(011^2,2(0) + Ci,3|(Al,3(0) , 



(29) 



must solve 



i-r\m) = {(Hl + Hr) + H,mt)) ■ (30) 
at 

We assume the difference between the time-derivative of 
l'A2.2(0) and (Hl + -Hfi)|(A2.2(0) is small (which assumes there 
is only a small amount of relative excitation), and neglect 
the time-derivative of |i/'3,i(0); by our initial assumptions, the 
prefactor C3j(0 is small. The time-derivatives of {c„,ai_„(0) 
are thus given by: 

;[C2,2(0|iA2,2(0) + C3,l(0l<A3,l(0>] - HiciMMt))- (31) 
Hence, [using that \^T,^\{t)) and |iA2,2(0) are orthogonal] 

iC2,2(0 - C2,2(0<«A2.2(0I^/I'A2,2(0), (32) 
iC3,l(0 - C2,2(0<iA3,l(0l^/l'A2,2(0)- (33) 

Within first order perturbation theory, {i/'2,2(0l^/|i/'2,2(0) is pe- 
riodic with a periodicity {T - n) half that of the oscillator 
period. The matrix element (i/'3,i(0l^/l'A2,2(0) is a product 
of a function with period T - n, and the complex exponen- 
tial exp(-/A£'intO of the energy difference between the intra- 
cluster degrees of freedom in both configurations (as plotted 
in Fig. 1). 

Denoting the periodic component of the interaction terms 
{4'n,N-n{t)\Hi\\l/2.2it)) as fn,N-n{t), wc must therefore solve 

/C2,2(0 ^ C2,2(0^/2,2(0 , (34) 
!C3,l(0 - C2,2(0^/3,l(OeXp(-!A£'intO , (35) 



with the boundary condition C2,2(0) - 1 . We first assume that 
the initial separation xq, and the interaction strength \g\, are not 
large. Within this regime we assume we can approximate fit) 
by a first order Fourier series f{t) a; 1 - cos(20, which implies 
all f„^N-n{t) differ only by a constant value; hence f„,N-n{t) - 
Af(t), /i,3(0 = Bf(t), with A and B dependent, in principle on 
g, and quite heavily on xq. We can use this to solve Eq. (35): 



C2,2(0 -exp / 



I dt' 

Jo 



Agfit') 



(36) 



; exp(iAg[t - sin(20/2 + . . .]), 



and if we neglect AEim under the assumption that the relative 
energy on both sides is similar. 



C3,l(0 ^ 



B 



exp 



i^igA dt' 



/(O -1 



(37) 



For short times, we can expand C3i(0 ~ B[igt + 
Oig^t^, g cos{2t))], i.e., proportional to gt and oscillatory 
terms and hence giving a linear increase when f - rni. At 
longer times the phase evolution of C2,2(0 becomes important, 
leading to cancellation in the terms of C3,i(0 and giving os- 
cillatory behavior with a period dependent on g. The linear 
increase with g after a collision is not expected to continue 
when § ^ 1 as higher-order terms become increasingly im- 
portant and the perturbation theory breaks down. 

We have so far neglected the difference in internal energy. 
This will introduce an additional phase between C3 1(0 and 
C2,2(0- With this included we have 



C3,i(0^ r ^/f'[/gBc2,2(f')/(f')exp(-/A£intf')] 
Jo 



^igB 



(38) 



) exp(/[Ag - A^int?' - osc]) , 



with "osc" denoting oscillatory terms such as k cos(2f), which 
are periodic with t — » t + n ox shorter fractions of n for 
the higher-order terms. Summing together terms of differ- 
ent phases will produce cancellation, hence if the exp(/[Ag - 
Aiiintf]) term has the same periodicity as f{t) and the "osc" 
terms, both n, the overall increase will be linear in time with 
no higher-order polynomial terms. This could therefore lead 
to resonant (suppressed) transfer if Ag - hEmt ~ n with n even 
(odd), and slightly suppressed transfer if n is a rational num- 
ber not close to an even integer, e.g. 1/2, 1 /3, 3/2. As noted 
earlier, the g — > 00 limit gives A/Sint ~ 1 and thus should lead 
to suppressed transfer if \Ag\ < 1 /2. We note that when \g\ ~ 
this resonance condition appears to be matched up to a factor 
g[A - (27r)"^], giving very long cancellation periods, however, 
as we see on Fig. 4 (and by the fact the perturbation strength 
scales oc g) the rate of atom transfer scales proportional to g 
and so cancellation can still occur before a significant popula- 
tion transfer is achieved. 

This simple analysis neglects higher-order effects such as 
pseudo-periodicity, and intra-cluster excited states are not 
treated explicitly. However, qualitatively we expect an ini- 
tially weak linear increase with long time oscillation effects 
for small \g\, and for ^ ^ 1 the timescale of these oscillations 
should drop. 
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D. Amplitude bound to oscillations 

One can look at each left/right number eigenstate [Eq. (19)] 
separately, assuming we have a probability of p for |2, 2), and 
of (1 - p)l2 for |3, 1) (with the same for the |1,3) state), no 
occupation of |4, 0) or |0, 4), and that there is no overlap be- 
tween the states and no mixing via the Hamiltonian. We can 
then state the energy £13 = (1, 3|//|1, 3) as follows 

£1,3 - £'pot,l + £'pot,3 + £'kin.l + £'km,3 + ^int.S ■ (39) 

Each term in this equation refers to the kinetic, potential and 
interaction energy of each side, with one or three atoms, re- 
spectively (note there is no interaction energy for the single 
atom side, taken without loss of generality as being left). Not- 
ing that the kinetic and potential energy terms must be posi- 
tive, we can derive the inequality 

^pot,! < -£1,3 - (£^11,3 + £'int,3) ■ (40) 

Using the conservation of E - {H) and t^E^ - {H^) - E^, it 
can be shown that (see Appendix D ) 



\Ei,,-E\< ^j^^AE , (41) 

which is equivalent to 

E-^^AE^E,,^E.^^AE. (42) 

Combining the upper bound of the above equation with 
Eq. (40), we obtain 



£pot,l <\E+ y ^^AfiJ - (£kin,3 + £mt,3) ■ (43) 

Finally, noting that Ej,ot,i - > <x>^/2, with the <(9>i 

meaning the expectation value of the 1 particle side of the 
wavefunction, we can obtain an inequality for the 1 atom po- 
sition expectation value 



<x>i < V2 ^iyE + ^jl^^E^ - ^mt.3 ■ (44) 

We can see that larger, positive g will constrain this bound, up 
to a point of saturation at the Tonks-gas limit, whereas poten- 
tially it is unbounded as g — » -oo (energies in this regime scale 
proportional to -g^ [29]) as the atoms gain a large amount of 
energy. 

V. POSSIBLE EXPERIMENTAL REALIZATION OF THE 
FOUR ATOM SYSTEM 

A. Optical lattice scheme 

Our results could be tested by creating an optical super- 
lattice [43], of two overlapping lattices, with one double the 




Displacement (a.u.) 



FIG. 2. (Color online): Potential from an optical super-lattice cre- 
ated by overlapping two lattices (all units are arbitrary). Grey circles 
represent a loading of two atoms in the ground state of each well. 
Our suggested scheme tunes the interactions to the desired value and 
then turns off the double frequency (dotted line) lattice leaving only 
the broader lattice (dot-dashed line), after which the atomic dimers 
collide. 



frequency of the other, then loading this with two atoms per 
site (in the ground state) in a Mott insulator regime [44]. This 
is shown schematically in Fig. 2. The interactions could then 
be tuned to be attractive via a magnetic Feshbach resonance, 
at such a rate that tunneling between sites is small, but the two 
atoms on each site tend to the ground state given by Eq. (11). 
The double-frequency lattice could then be ramped down, 
leaving only the wider lattice, thus creating the initial con- 
ditions of two equally separated dimers in an approximately 
harmonic potential. 

Some freedom with xq could be achieved by modifying 
the relative strength of the double-frequency lattice compared 
with the primary lattice. Reducing it will push the minimum 
closer together, but also make tunneling between the sites 
more significant. Careful ramping-down schemes of the laser 
power of the double-frequency lattice could also be incorpo- 
rated, which would give further freedom to move the sites 
closer together after creating the dimers. Slower ramping will 
also make things closer to adiabatic, thus reducing the excita- 
tion in each dimer created by the switch-off. The relative ve- 
locity between the two dimers in terms of the final harmonic 
oscillator units will equate to an effective initial separation — 
approximately the separation the dimers will reach after the 
first collision. A faster (slower) ramping scheme would give 
a larger (smaller) effective xq, however, to be most applicable 
with the results of this paper, a slow scheme would be ideal to 
minimize excitations and minimize the degree of anharmonic- 
ity in the potential that the dimers sample. 

After some free-evolution time, the double frequency lat- 
tice could then be quickly restored with an extremely high 
lattice depth, separating the left and right components of the 
wavefunction, with no further tunneling possible. This would 
allow for a direct measurement of as defined in Eq. (13), by 
then imaging the lattice with resonant light; light-induced col- 
lisions [45] will reduce this to a parity measurement with an 
empty site being either a zero or two population, and a single 
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atom being a one or three population. This is actually suffi- 
cient information, assuming we know the total atom number 
in the two sites was exactly 4. In terms of the states given in 
Eq. (19): no atoms on either site is a measurement of a |2, 2) 
configuration (or a |4, 0)/|0, 4) configuration, but this is only 
significant during collisions), both sites occupied is a mea- 
surement ofa|3,l)/|l,3) configuration, a single occupied site 
and an empty site would imply some inelastic process has oc- 
curred (such as three-body recombination or background gas 
collisions) and such a result would thus be null. 

If the effective xq were an appreciable fraction of the lat- 
tice width, this scheme could also show some more interest- 
ing physics beyond the scope of this paper, with collisions 
coupling energy into the center-of-mass mode and the tunnel- 
ing of the single atom in the single-trimer states (considered 
in Sec. IV C) to adjacent lattice sites. It could even have a ki- 
netic energy greater than the maximum barrier height between 
sites and join an effective conduction band [46], allowing for 
entanglement between lattice sites. These efifects may also be 
worthy of experimental investigation. 



VI. NUMERICAL METHOD 
A. Basis set expansion 

To perform many-body computations we expand the field 
operator over the set of Hermite functions of a given width W 



<Pk(Wx) = y HkiWx) exp , (46) 

with Hii(x) the Hermite polynomials, and diagonalize the 
Hamiltonian in a Fock state basis \no, . . . ,«t)o), truncated via 
the condition 2^ kn^ < rj. Such a calculation would require 
an unfeasible amount of states to converge, were it not for the 
fact that the center-of-mass part of the Hamiltonian commutes 
with the rest of it. This means we can just consider a subset 
of this truncated Fock space where the center of mass of the 
gas is in the same state. This does not have to be the ground 
state, as we simply ignore the center-of-mass time evolution 
and can account for it later. The procedure essentially involves 
diagonalizing the finite basis in terms of the operator 



B. Experimental parameters 

In terms of typical experimental parameters, the .s-wave 
scattering lengths would need to be very substantial in order to 
give measurable effects. Strong interactions generally require 
tuning scattering lengths near to Feshbach resonances, and in 
such strongly interacting regimes confinement effects can shift 
the effective ID scattering length if aja^ is not small [32]. 
The chosen Feshbach resonance would ideally be broad, min- 
imizing uncertainty in the effective interaction associated with 
a lack of precise control of magnetic field fluctuations. 

Alternatively, some atoms such as Cesium can have large 
"background" scattering lengths far from resonances [47], 
e.g., ~ +3000ao where uq =s 5.3 x 10""m is the Bohr 
radius. In terms of a rescaled g parameter in harmonic oscil- 
lator units, if one had ~ In x IHz and very strong radial 
confinement a)_i_ ~ 2nx 0.4kHz, we have 



g^loj^a, ±1.2, (45) 



which is of unitary order. 

We essentially have three experimentally tunable parame- 
ters, Us, (±>x and (jj^ which can be varied smoothly with small 
adjustments to a magnetic field or modifying laser powers , 
focusing, or detunings. However, dropping is undesirable 
as it increases experimental timescales, and increases the like- 
lihood of background gas collisions; additionally, unwanted 
three-body recombination effects scale oc \ast (generally be- 
ing worse for < 0) meaning one would need to determine 
an appropriate compromise solution. 



A"'"A = 2 V(^ + 1)(; + l)fll+iflAfl]flj+i , (47) 
K) 

(where A' — 2^. + la^^jat is the creation operator for a 
dipole mode of width W) and taking the eigenvectors with 
eigenvalue zero; this procedure is discussed in detail in [28]. 
For the calculations in this paper we use the eigenstate width 
W = 1 as the harmonic oscillator length is always a relevant 
scale. 



B. Convergence testing 

We first need to represent our initial condition in terms of 
this basis set, noting that due to the truncation the state cannot 
be represented exactly, with larger initial displacements and 
larger interaction strengths harder to represent in this basis. 
We require a reasonable fidelity of our numerical initial con- 
dition to the true state, achieving fidelities of > 99.5% for all 
the numerics used in this paper. 

Measuring convergence during time evolution with such a 
method is more difficult. Performing the calculations with a 
variety of basis sizes and calculating the fidelity over time can 
give an indication for how long the calculations are reliable, 
for which we plot, in Fig. 3, our most extreme values of g. 
This is probably the strictest measure of convergence applica- 
ble, given the large number of degrees of freedom in a many 
body wavefunction, for example a product state with a large 
number of atoms would have a fidelity exponentially tending 
to zero for any finite difference in the product wavefunction. 
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FIG. 3. (Color online): Plot of \{t//(v, 0|iA(Tmax, 0>P; the fidelity of the 
wavefunction computed with smaller basis (energy cut ofi' at v) to the 
wavefunction computed using a larger basis truncated at Vn,„ = 113. 
This shows the extreme values of g employed in the numerics, lower 
absolute values of g converge more rapidly. 



VII. NUMERICAL RESULTS 
A. Preamble 

All the results graphed here are calculated for N - A and 
jcq = 3 in order to investigate the elfects of varying interac- 
tion strength for small numbers. In general smaller xo greatly 
increases interaction times between clusters and thus rates of 
atom transfer It also reduces the amount of free energy in the 
system, however a greater amount of the wavefunction will 
be found towards the center at all times and thus expectation 
values of .^Vr will be harder to interpret. The results here are 
broken down into three sections, the first examines the vari- 
ation in left right number, the second examines the variance 
is position about one side and the final section examines the 
single body von Neumann entropy. 

B. Left and right particle number dynamics 

Because our initial condition has a definite number of two 
atoms either side of the trap, the left-right number uncertainty, 
AA^«, in our system is initially very near zero. We note that a 
mean-field-like state or a symmetric superposition of 3 and 1 
atoms either side both give AA^r = 1 , which is also the value 
this quantity will take in our non-interacting system when 
each of the clusters collide. We therefore first consider the 
minimum to minimum values taken by ANr before and after 
each collision. The change after the first collision is given in 
Fig. 4 and the change over the first 150 collisions is plotted 
in Fig. 5. Despite the fact that the increase after the first col- 
lision is similar for both attractive and repulsive interactions 
of similar magnitude, the long time change is very dilferent, 
with the timescales being much longer in the attractive case. 

In either case, the left-right number does not reach an equi- 
librium on the timescales considered, with oscillations and 





Interaction strength: g 



FIG. 4. (Color online): (a) Minimum value taken by ^.Nr [Eq. (B6)], 
after one collision (b) Frequency difference of peaks in the Fourier 
transform of KNr from the non-interacting values {t = nix) divided 
by n. (a) shows that for g > Q, the increase to number uncertainty 
is greatest for g x 2.3 and decreases when interaction strength is in- 
creased further. The g < behavior is initially similar but deviates 
at around |g| = 0.6; rather than saturating it appears to increase even 
more rapidly with |g|. It is not clear what will happen for g < 
and Igl » 1, which will be a topic for further investigation, (b) 
demonstrates the existence of pseudo-periodicity in the system (in 
addition to low frequency components relating to the long time be- 
havior). The non-interacting system has frequency peaks at f„ = n/n, 
the quadratic fit (solid line) indicates these peaks shift by an amount 
roughly equal to -ng/lOOn. 



revivals present. The time-dependent perturbation theory of 
Section IV C indicates that atom transfer processes are sup- 
pressed by an internal energy difference between the |2, 2) 
and |3, 1) configurations of the wavefunction, which leads 
to destructive mixing over a few collisions, unless a phase 
matching condition occurs. If intra-cluster excited states (dis- 
cussed in Section IV B 3) are present, the energy difference 
between each configuration, A/iiei, may be small (along with 
Ag) meaning cancellation occurs on longer timescales, lead- 
ing to fluctuations in ANr over 10s of harmonic oscillator pe- 
riods. 

Figures 6 and 7 (a) show the amplitude of each number 
component in the wavefunction as it evolves in time for g = 3 
and g = -1.7; note Fig. 5 takes only the minimum values of 
these curves to avoid the spikes on collisions. The maximum 
amplitude of the |3, 1) and |4, 0) components (at least initially) 
occurs on collisions (corresponding to a minimum amplitude 
of |2,2)). Decreasing of this peak amplitude may be inter- 
preted as the time of collisions between clusters becoming less 
well defined, due to the distance between their centers of mass 
becoming less well-defined (i.e., its coiTesponding probability 
density becomes broader) and the forming of intra-cluster ex- 
citations. 

At late times (f > 100) on figure 7, all the expectation val- 
ues for n 7^ 2 are almost the same as those for Gaussians cen- 
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tered on zero. This is due to only the two-dimer (attractive 
n - 2 ground states) setup being significant, as the exciting 
of "intra-cluster" excitations is suppressed by the large energy 
gap, and atom transfer interactions are suppressed by an en- 
ergy difference, leading to a phase mismatch and hence a can- 
cellation. However, energy is still transferred to the relative 
position wavefunction (described in Section IV B 2), increas- 
ing the uncertainty in the separation of dimers, and so some 
component of the wavefunction is always undergoing a col- 
lision yielding a finite value for the left-right number uncer- 
tainty. As a result of our scaling in Eq. (17), the n 2 values 
are just those of the dimer system in collision, and only a small 
contribution to |3, 1) comes from states that are similar to a su- 
perposition of a cluster of 3 atoms to the left (right) and a free 
atom to the right (left). 
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FIG. 5. (Color online): Minimum value obtained by ^Nr, as given 
by Eq. (B6), after a given collision. For weak interactions (|g| < 
0.1) the behavior is the same for attractive and repulsive, but for 
slightly larger values there is a clear difference in the timescales, with 
repulsive interactions producing larger number uncertainties more 
quickly, despite the fact that Fig. 4 shows there is little difference in 
AA'k after one collision. This difference is likely due to the increased 
(decreased) energy spacing between the ground and first excited state 
of the two atom system with attractive (repulsive) interactions, dis- 
cussed in Section IV B 3, and the energy difference between the two- 
two and three-one number configurations, as discussed in Section 
IV C, which leads to a phase mismatch. For large repulsive values 
(g > 2), AA'r reaches a maximum value and then undergoes complex 
partial revivals on timescales of 30 time units, (tens of collisions). 



C. Equilibration of energy into inter/intra-cluster excited 
states 

We wish to quantify the amount of energy transfeiTed from 
the center-of-mass energy of each cluster to excitations be- 
tween the atoms, as discussed in Section IV B 2 and Section 
IV B 3. We therefore investigate the standard deviation in the 



FIG. 6. (Color online): For g = 3, xo = 3, (a) time evolution of 
probability of finding n(oT N - n) atoms to the right with the ampli- 
tudes of wavefunction components decomposed into eigenfunctions 
of number L/R number operator, (b) Expectation value of position 
to the right on sections of wavefunction decomposed into eigenfunc- 
tions of L/R number operator, (c) Variance in position to the right 
as defined in Eq. (48), paralleling (b). The expectation value to the 
right [(b)] effectively tracks the particle-like motion, but after long 
times the motion appears effectively damped, (c) can quantify this 
effect — the peaks of tr(n, N-n) increase from their initial value and 
continue to oscillate about a maximum, except for o"(4, 0) (which is 
only significantly probable during collisions) indicating a transfer of 
energy to the degrees of freedom described in Sections IV B 2 and 
IV B 3. This remains true even at very long times t ~ 1000, with 
progressively smaller partial revivals and so can be said to have equi- 
librated. 



position to the right, for a given number of atoms to the light 

o-(n,N - n, t) - 

^mxW^{x)x^^{x))f - \{®{xn\{x)xMx))f\'^ , (48) 

essentially the width of the atomic density distribution on the 
right hand side, about the expected value for position, given 
that n atoms are on the right-hand side. This is plotted in 
Fig. 6 (c). The repulsive case shows a consistent increase in 
the height of the peaks (excepting the n - A peak), with only 
small periodic oscillations. The attractive case however shows 
cr(n, N-n,t) to be initially similar but then dropping to a min- 
imum value for n + 2. We note cr(n, N -n,t) cuts off anything 
on the left side, and so is difficult to relate to the amount of ex- 
citation if the left and right states are separated by a distance 
smaller than the size of their internal structure, as they will 
contribute to all the n + 2 expectation values. Intra-cluster ex- 
citations as we have defined them are present if the wavefunc- 
tion either side of the center does not look like a displaced n 
atom ground state; it is possible such excitations could reduce 
the position uncertainty but are generally expected to make it 
broader and thus increase cr(n, N-n, t). These excitations are 
dominant processes in the increasing of cr for the repulsive 
case plotted in Fig. 6 (c), and appear to persist at long times. 
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FIG. 7. (Color online): The same quantities as Fig. 6 but with 
g = -1.7. The short-term behavior of the expected one-atom po- 
sition (b) is similar to the repulsive case, but is increased in magni- 
tude. At long times, the right-position expectation values drop to an 
approximately constant value for all but n = 2, this being the value 
of a Gaussian state in the center of the trap, for reasons explained 
in section VII B. This is also the case in (c) — essentially the only 
significant contribution to the n t 2 states comes from uncertainty 
in the separation of the atomic dimers, which smooths over transfer 
effects. 



For the g < case, at very early times, say t < 20, the 
contribution to cr(3, 1, t) from states in the single particle and 
cluster-of-3 configuration is visible. By (approximate) mo- 
mentum conservation the single atom must have considerably 
more energy after a collision than the 3-atom state, which ex- 
plains the large n - I position expectation values away from 
collision. However, in the strongly attractive case this transfer 
process is cycUc, and it never transfers large populations to 
these configurations. As we noted before, contributions can 
come from an oscillating dimer state if the relative separa- 
tion is small. Initially this only occurs during collision, but 
inter-cluster excitations (which can be interpreted as an in- 
creased uncertainty in how much the centers of each cluster 
have shifted due to interactions), lead to an increase in rela- 
tive position uncertainty^. Hence, at late times there is always 
significant wavefunction density in the trap center, that is to 
say at any time t > fiate some non-negligible pait of the wave- 
function is always undergoing collision. Hence, if the contri- 
bution from the singlet-triplet state is too small to see we can 
conclude that the cr(2,2,r) reaching a maximum corresponds 
to this mode reaching a steady configuration. This is the dom- 
inant effect in the attractive case shown in Fig. 7, but is also 
present for g > 0. 



^ although Fig. 5 indicates this process undergoes partial revivals 



D. Relaxation to equilibrium 

One questions of interest is whether the system reaches an 
equilibrium at long times. We attempt to quantify this by look- 
ing at the single body density matrix and its von Neumann en- 
tropy, given by Eq. (21); however, this quantity (like most in 
our system) has a time-dependence due to the repeated colli- 
sions that are a consequence of the system as a whole being 
held within a harmonic confining potential. In order to sim- 
plify our analysis we look at the time averaged value over a 
period of T - In and quantify the degree of short-time change 
via the variance of this average. These are plotted in Fig. 8; 
(a) shows that for both positive and negative g, S vn increase 
towards a maximum value, with small amplitude oscillations 
in a similar way to AA^ but with much smaller variations. For 
fixed \g\, the ^ > entropy generally increases slightly faster 
and to higher values than the equivalent g < case, but is oth- 
erwise quite similar Fig. b) shows the standard deviation over 
the 27r averaging period, the rapidly changing (time scales of 
less than 27r) effects continue for much longer in the attractive 
case compared to the repulsive. Transfer effects [discussed in 
Sec. IV B 4] are likely the cause of this short time oscillation 
as they are predicted to be cyclic on the timescale of a few col- 
lisions when g ~ 1. The variation dying down at long times 
can be explained for the ^ > case by intra-cluster exited 
states breaking the cyclic effect, and for ^ < 0, by the slower 
effect of the broadening of the inter-cluster wavefunction to 
the point where the collision time is not well defined. 



VIII. CONCLUSIONS 

We have considered a system of A = 4 atoms with contact 
interactions, confined within in a harmonic potential. Our ini- 
tial condition was a symmetric setup of two A/2 atom ground 
states, displaced from one another by a distance jcq (taken to 
be 3 harmonic oscillator lengths for most of the numerics), 
which we then left to oscillate and undergo collisions. Ini- 
tially there is no entanglement between the atoms on the left 
and on the right, however interactions lead to the generation 
of entanglement. 

We investigated left/right number variation within the sys- 
tem, based on an operator which could in principle be mea- 
sured directly in the experimental setup we suggest in this 
paper Initially both (left and light) states have a near def- 
inite number of two atoms and hence a number uncertainty 
AAr, which is initially close to zero. When the left and right 
states are well separated, AAr is a measure of entanglement 
between the left and right sides. However when the two states 
are close, i.e., during collisions, AAr ~ A/4 = 1; we therefore 
investigated the difference from minimum-to-minimum value 
taken over a time range of around n, i.e., the minimum value 
of AAr obtained after the nth collision. There is a marked dif- 
ference in the evolution of AAr between the g <0 (attractive) 
and ^ > (repulsive) cases. When \g\ ^ 0.5, number uncer- 
tainty builds up much more slowly with attractive interactions 
than with repulsive, essentially resisting entanglement. This 
is despite a large increase to the change in number uncertainty 



13 




50 100 150 200 250 300 




50 100 150 200 250 300 

Time (harmonic oscillator units) 



FIG. 8. (Color online): For jcq = 3, (a) von Neumann entropy aver- 
aged over a time period of 2;r, as defined by Eq. (21) and Eq. (22) 
(b) the standard deviation of this quantity, given by the square-root 
of Eq. (23), for a range of interaction strengths both repulsive and 
attractive. Entropy increases gradually at early times / < lOn, then 
increases at a more rapid rate before leveling off to an almost con- 
stant value with small fluctuations. This behavior is similar for both 
attractive and repulsive interactions. The variance over the 2n aver- 
aging range behaves very differently for strong attractive and repul- 
sive interactions, with the short-timescale fluctuations persisting for 
much longer if g < 0. This difference is explained by a change in 
the dominant processes, with the attractive system being unable to 
excite the relative degrees of freedom in a cluster and thus transfer 
of atoms between each cluster becoming more significant. Fig. 6 (b) 
shows atom transfer dynamics in the repulsive case have only small 
fluctuations at late times. 



that is generated by a single collision. This increases quadrat- 
ically with \g\ when g ^ -1.3, but in the repulsive case the 
increase reaches a maximum, and then drops as g increases 
further Additionally for ^ > we observe long-timescale 
high-amplitude number fluctuations, which continue even at 
late times (over 100 collisions). 

This behavior is explained by our time dependent perturba- 
tion theory on the atom transfer process, and the energy differ- 
ence between the intra-cluster excited states. We investigated 
the effect of AEi„t, the energy difference in intra-cluster ener- 
gies between the {2, 2) (two displaced N - 2 ground states) 
and {3, 1) (one free atom and one N - 3 atom ground state) 
configurations. Assuming the average interaction energy be- 
tween the clusters to be weak (i.e. \Ag\ <sc 1), increases to 
lAEintl lead to a phase mismatch and thus to destructive in- 
terference so that the population transfer cycles periodically. 
If intra-cluster excited states are present, this picture breaks 
down, since each of these excited states phase-evolves at a dif- 
ferent rate; cancellation becomes more complicated and the 
states less localized, which occurs for large g > at long 
times. The energy gap between the ground and excited states 
of each of the Njl atom clusters is increased (decreased) when 
g gets smaller (larger), which reduces the maximum popula- 
tion that can be transferred to excited states. The excited states 
become effectively inaccessible as g <k 0, resulting in an ef- 



fectively two-level system of the {2,2) and {3, 1) configura- 
tions. Our perturbation theory indicates that for sufficiently 
strong attractive interactions, with very specific values, phase 
matching would be possible, allowing for resonant transfer 
However this is outside the regime our numerical method is 
capable of reliably portraying, and will remain an avenue for 
future research. 

By separating the system into components of the wave- 
function with definite number (number states of the number- 
to-the-right operator) we have observed the evolution of the 
positions associated with one/two/three atom number states, 
and the right side position variance. For § = 3 the peaks in 
position variance increase to a maximum for all Nr - n m 
around 100 harmonic time units (100/2;? oscillator periods or 
around 30 collisions) and do not fluctuate greatly. Consider- 
ing instead the case where g - -1.7, after 60 collisions, we 
find that for Nr 2 position and position uncertainty are the 
same as they are for a state undergoing collision, whereas the 
Nr - 2 tends to a maximum. This indicates that the state is 
well described by two atomic dimers with a significant un- 
certainty in their relative displacement and almost no ampli- 
tude of a singlet-trimer like state is present in the wavefunc- 
tion; this motion again undergoes partial revivals on very long 
timescales. 

In addition, we have investigated the von Neumann entropy 
of the single-body -density matrix 5vn(0> in order to inves- 
tigate to what degree the system tends to an equilibrium. We 
note 5 vn(0 is zero for a product state (all atoms with the same 
wavefunction/occupying the same mode) and can be consid- 
ered a measure of how mean-field-Iike the state is. Addition- 
ally 5vn(0 is constant for our system if g - Q, despite the 
wavefunction evolving periodically in time. At long times 
with repulsive interactions, 5 vn (time averaged over a period 
of 27r) increases to a steady value with only small fluctua- 
tions over the averaging period. However, long-term fluctu- 
ations (over the order of twenty k time units) are still present 
and appear to be due to atom transfer processes which do not 
appear to equilibrate on the timescales considered in this pa- 
per The time required to reach maximum entropy decreases 
with larger g but this appears to saturate with little change for 
g ^ 2; for an initial separation of xq - 3 this takes around 
30 collisions. This short-term increase appears to be due to 
the inter-cluster degrees of freedom discussed in the previous 
paragraph; the associated probability density with the separa- 
tion of the two clusters becomes less peaked. With very weak 
attractive interactions, the system's behavior is similar to the 
repulsive case, however for \g\ ^ 0.5 higher intra-cluster ex- 
cited states become less accessible, leading effectively to a re- 
duction in the number of accessible degrees of freedom, such 
that the left/right states behave more like solitons. In this case, 
the time average of 5'vn(0 does not tend to a long-term mean 
value as compared with the case of repulsive interactions of 
similar magnitude; there is also a great deal more short-time 
variation, which persists for longer. The shoit-time variation 
can be attributed to the strong atom transfer effects, which are 
predicted to cycle population continually due to an energy dif- 
ference. The effect eventually reduces as displacement uncer- 
tainty between the two bound states (which now behave like 
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quantum solitons) increases, which is the mechanism behind 
the long term entropy increase. 

A pseudo-periodicity effect is also present. The non- 
interacting system is periodic with a period tt, and thus the 
Fourier transform of any time dependent expectation values 
will have frequency peaks at n/n. We have examined how 
these peaks shift for the left/right number uncertainty as inter- 
action strength is varied and have found an approximately lin- 
ear shift with g over the range considered. Changes to higher 
order components of the frequency spectrum depend deviate 
slightly from the linear dependence shown by the first order, 
with differences only clearly manifest for \g\ > 1). 
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Appendix A: Identities involving Jacobi coordinates 

1. First identity 

We wish to show that the Jacobi coordinates defined by 
Eq. (5) and Eq. (6) satisfy 



2. Second identity 

We rephrase Eq. (6) as = + - 1)1 2*= j -^^j- 

Recursively substituting in equivalent expressions for 
Xk-i,Xk-2i ■ • ■ > xn/2+\ yields (for N/2 + 1 < k < N) 

k-\ ^ J Nil 

Z 7 + ]w22"- ^^^^ 

;=/v/2+i J '^'^ j=i 

and for k - N/2 -H 1 we have Xff/2+1 - ^n/2+\ + (2/N) Tj^li ^j- 
Hence, summing over all k e {N/2 + l,N/2 + 2,..., N], 



N k-l ^ N/2 



Z Z z z ^z 

k=N/2+l k=N/2+l k=N 12+2 i=N l2+\ ■' k=\ 



Xk 



N-k 



k=N /2+1 k=N /2+1 k=l 

N N/2 



Z 



Xk 



k=N/2+\ k=\ 

from which we deduce the desired identity: 

N N/2 N 



k=N/2+\ k=l k=N/2+\ 



(A7) 



Z-^ 



^^C(N) + 



k=\ 



z 

k=2 



k-l . 



(Al) 



We prove this inductively. The N - 2 case can readily be 
verified, after which we may consider the increase of number 
from - 1 to A^. In particular. 



N-\ 



Z 4 = 4 + - 1)4(^-1) + Z 

k=\ k=2 

Noting that - x^ - xc(a'-i), we then deduce 

N 

Y^xl + (A? - 
N-l 



k=l 



-,2 A: — 1 . 

[xn - xc(iv-i)J + 2^ —j—^i 



N 



k=2 



Collecting terms, this reduces to 

N ^ N ^ 



k=\ 



k=2 



=NXc^f,y + 



z 

k=2 



k- 1 , 



(A2) 



(A3) 



(A4) 



which completes the proof An equivalent result also holds in 
3D [48]. 



Appendix B: Calculations for the number- to-the-right operator 
1. Analytically determined properties of 

From the definition of Eq. (13), it follows that 

A^^ = I dxdx^>\x)¥(x)'i'(x)^'(x) + NR, (Bl) 
Jo 

and, given a general (symmetrized) many-body wavefunction 
i/>{x), one may deduce the expectation values 



r) ^N I dxi I dx2... dxN \(f/(x)f- , (B2) 

Jo J-co 



(Nl) ^N(N 



- 1) I dx\dx2 I 
Jo J-i 



dx\dx2 I dxT, . . . dxN \i/f(x)\'^ + (Nr) ■ 



(B3) 



For a product-state wavefunction (f/{x) = nf=i '/'(^k), expecta- 
tion values are simple to calculate, as all integrals are separa- 
ble and most evaluate to unity. In this case 



r}^N [ 
Jo 



dx |0(x)r , 



{Nl}^N(N-l) 



f 

Jo 



dx \(p(x)f 



+ {Nr} 



(B4) 



(B5) 



= [{N-l)/N]{NRy- + {NR}, 
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and so the variance of for a product state simplifies to 

ApNr = <^R>(1 - {Nr)/N) . (B6) 

We may determine analytic expressions when g = 0, which, 
for the purpose of this paper, we limit to the N - A case. 
Without interactions, our many body wavefunction is given 
by Eq. (9), and 

dx\^{x, t)f- = T [1 + erf(ji;o cos(f))] , (B7) 
2 



3. Numerical calculation of restricted region expectation 
values 

In addition to this we wish to calculate expectation val- 
ues in restricted regions via Eq. (17), corresponding to sec- 
tions of the wavefunction with exactly n particles to the left 
or right, along with the associated normalization factors when 
the wavefunction is divided into these regions. If our many 
body wavefunction is (/'(x) then the normalization factors are 
given by 



^..0*(.,±xo,O0(x,..o,O = e-o-«-<-)^ (B8) ^„.^^ r..,.....„ f ........... 



n * 1 

I dxd) (x, ±xo, t)4>(x, +xo, f) =r- [1 + erf(;ico sin(?))] 
Jo 2 (B9) 

^ g-jr5±ixo sin(20/2 

with erf denoting the error function. Calculating (.^V^) in prin- 
ciple requires accounting for 36 different terms, however, as- 
suming we can neglect terms proportional to exp(-2xQ), only 
6 are important, and we have 



{[l-erf(xocos(f))]' 



24 

-i-4[l - erf^(xQ cos(f))] 
+ [l+erf(xocos(f))]'} + <^R) 
-5 - erf^[xocos(f)] - 



(BIO) 



Subtracting 4 then yields the variance as given by Eq. (16). 



2. Numerical calculation of number variance 



(B13) 

and the expectation value of the distance to the right operator 
is equal to 



{x 



dxi... dxN ^ Xk0(xk) 

k=o 



--N: 



Y]&(x,) Y] @{-xj)\ii^{^f 

P k=l j=n+l (B14) 

n)!w! Jo 



{N- 



Ml N 

: I dx„^i...dxN ^ xi,\if/(x)\^ . 



For computation, these operators are converted into matrix 
form by taking the matrix elements between different ele- 
ments of the basis set, and then projected to our reduced 
(center-of-mass ground state) basis. 



In order to calculate the number variance we decompose 
the field operator into our basis set, ^(x) - ^k<Pkix)- In this 
form we can express .^V^ as 



A'r = ^ yikyj(a\a\akai + Nr , 
ij,k.e 



(BID 



where yjc — dx (pj(x)(pf(x) is the positive space overlap 
between two Hermite functions, given by dj(/2 if 7 + ^ is even, 
and otherwise given by 

yje = (-l)(>+^-l>/2 J _ _ ^/2; 1 - + f]/2, -1) 

2-j(j + { + 2)\\ 



(B12) 



where 2^1 denotes a standard hypergeometric function. 
Likewise the integral from minus infinity to zero is (-ly^^yjc. 
This formula is useful for small numbers and testing, but for 
practical purposes we calculate the integral via Gauss La- 
guerre quadrature, which is numerically exact for odd j + £ 
(all other cases are trivially zero or one half) given a rule of 
order (j + { + l)/2 or higher. Given our truncated basis and 
symmetry about x = 0, this can be expressed as a finite size 
matrix of only even-parity functions with {Nr} - A^/2 just a 
numerical constant for our initial condition. 



Appendix C: Two cluster wavefunction evolution 

Here we derive the time dependent wavefunction describ- 
ing the center of masses of our two cluster system, i.e. the 
part acted on by Z/^^^, the center-of-mass components from 
Eq. (24); with H; ignored. Denoting yuyi as the coordinates 
of the center-of-masses of each cluster, up to a normalization 
factor our initial two-cluster wavefunction is given by 

n2\ 



<3'i.}'2l'^«,A'-«(0)> ocexp 



N- 



>'2 + 



N-n 



xexp(-^[yi -X„]2j-Hperm , (CI) 

with "perm" indicating the term obtained by permuting yi 
and y2, as required by symmetry. This gives rise to a time- 
dependent normalization constant which we do not discuss 
here. If we instead express this in terms of yc = [nyi + {N - 
n)y2]/N and = yi - y2 we have 

{yC,yR\<Pn,N-n(0)} CX exp ( 



2N[N - n] 



X exp 



Nyf^ 



+ perm , 



(C2) 
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where in this case "perm" is simply flipping the sign of jr, and 
we can factor out the yc dependence. If we temporarily ignore 
interactions between the two clusters, it is straightforward to 
generalize this to the time dependent case via Eq. (10): 



<3'c,3'rI'^«,a'-«(0) 0^ 

/ n[(N-n)yR-NX„cos(t)f 
exp exp 



( Nyl) 

I PYD 

2N[N - n] 

X exp (/ t - ny^X„ sin(f) H — -in ) sin(2f) ) 

\L 4 \ N - nJ IJ 

+ perm . (C3) 

Interactions between clusters can modify only the depen- 
dent part of this wavefunction. 

Appendix D: Energy bound for Hamiltonian variance 



As the Hamiltonian is time independent, the time evolution 
operator commutes with all powers of the Hamiltonian. De- 
noting our state as \tfr{t)} we have for any time t 



m)\H"m)) = {m\H"\m) , « = i, 2, . 



(Dl) 



As absolute values of energy are not physically important, we 
consider a re-zeroed Hamiltonian 



9<^H- {m\H\m) 



(D2) 



as it will make the mathematics more convenient. Introducing 
the notation for the variance of the re-zeroed Hamiltonian 



AE- = (-K^) 



(D3) 



we note that this quantity must be positive and real as is a 
Hermitian operator 

Let us define two wave functions \4'i{t)) and |i/'2(0) as being 
negligibly mixed at a certain point in time if 



(D4) 



with T] a small parameter. Note that in lattice models rf could 
be exactly zero up to some finite power n. If both the initial 
wave function and \4'\,i{t)) are normalized to one and the lat- 
ter are negligibly mixed, the wave function at time t can be 
written (up to a global phase factor) as 



m)) = VpI'Ai(o> + Vw^?''V2(o> (D5) 

with real a and < /:> < 1 . Introducing the notation 

{<ik")j = {^j{m"Wj(t)) , (D6) 

we can see from Eq. (D4) and the fact that the expectation 
value of total Hamiltonian is zero, that these two quantities 
are related via 



= —mi + o{ii) . 



(D7) 



Setting 77 = in Eq. (D4) we have for n - 2 



>pWl + (i-p)ml 



(D8) 



with the second step true again by the fact that H is Hermitian. 
Finally, substituting in for via Eq. (D7) we obtain 



AE"- > 



2 - ^ ^mi 



AE^ > 



i-p 



AE^->^ ' 



(D9) 
(DIO) 

(DID 



P(l-P) 

which leads to Eq. (41) in the main text. 

1. Analytic calculations of AE 

For our two particle initial condition, if ;ico » 1, i.e., well- 
separated initial clusters, we can analytically determine E and 
AE. Within this well-separated approximation we only need 
to consider one cluster, displaced a distance xq from the cen- 
ter, and multiply by 2 to get the values for the whole wave- 
function. For dimers, our wavefunction is f{xi - xo,X2- x^f^'' 
as defined in Eq. (11), otherwise it is not analytic. This wave- 
function is still an eigenstate of the relative Hamiltonian (for 
n particles), with some eigenvalue E^J^, but not of the center- 
of-mass part. Therefore we need only consider the center-of- 
mass Hamiltonian 



Hcixc) 



2n dx^ 



acting on the displaced ground state 

\l/4 



i/'c(JCc) 



expi-n[xc- xo]/2) 



(D12) 



(D13) 



to get all contributions to the variance. Acting the Hamilto- 
nian on this wavefunction we obtain 



HciJ/cixc) = 
Hl<Pc(xc) = 



v2\ 



• nxQX + 



nxQ 



lAc(JCc) 



- -I- — (4x - 3xo) + n^xlixQ - 2xf 



(D14) 

lAc(-^c) , 



which can then be used to determine the expectation values 

Ji 



1 «x„ 
{He) — , 

{Hl)^^-+nxl+"-^ 



(D15) 



AE can then be calculated as the standard deviation of two 
times He 

A£ = 2 .^{Hl) - {He? = V2^-Xo , (D16) 
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which is twice the square root of the difference between the 
initial (dimensionless) potential energy and the ground state 
energy. The reasons for this are similar to why a classical co- 
herent state with an average value of photons has a shot 
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